3d Interpolation Using Hankel Tensor Completion by Orthogonal Matching Pursuit

Introduction. Seismic data are often sparsely or irregularly sampled along one or more spatial axes. Irregular sampling can produce artifacts in seismic imaging results, thus multidimensional interpolation of seismic data is often a key processing step in exploration seismology. Many solution methods have appeared in the literature: for instance in Spitz (1991), Sacchi (2000) it was proposed to perform seismic trace interpolation that handles spatially aliased events using the linear predictors to estimate the missing traces through linear filters in the f-x(y) domain. Conversely, multidimensional Fourier reconstruction methods exploit a different data representation domain, while still assuming that seismic data consists of a superposition of plane waves: statistical sparsity is assumed in order to retrieve a model that consists of a few dominant wavenumbers representing the observations (Liu, 2004; Sacchi, 2010). In Xu (2004, 2005) and Sacchi (2000) the so-called Antileakage Fourier Transform (ALFT) algorithm was proposed: it tries to resolve, or at least attenuate, the spectral leakage of the irregular Fourier Transform by taking advantage of the compressive sensing framework (data are assumed to be sparse in the frequency-wavenumber domain). In more recent years, methods based on alternative local transformations, such as curvelet and seislet transforms (Herrmann, 2010; Fomel, 2010), have also been proposed. Recently, new interpolators have been developed recasting the interpolation problem to a compressive sensing matrix completion problem (Candes, 2009; Yang, 2013; Herrmann, 2014). A matrix with randomly missing entries can be completed by solving the rank minimization problem under the low-rank assumption of the underlying solution, that is in opposition to the fact that the subsampling operator tends to increase the rank of the matrix. Trickett (2010) proposed to apply a Cadzow filtering to solve the trace interpolation problem. The Cadzow algorithm replaces the block Hankel matrix of the incomplete and noisy multidimensional observations: it allows to recover the complete data volume by its low rank approximation and successive averaging along the anti-diagonals of the rank-reduced Hankel matrix. In recent literature, the low-rank tensor completion technique has become widely used in many areas such as computer vision, signal processing and seismic data analysis. This approach is based on a high order generalization of the low-rank matrix completion problem, and can be approximated by convex relaxation which minimizes the nuclear norm instead of the rank of the tensor. In Trickett (2013) a generalization of Hankel matrices is applied to the tensor to perform an efficient completion. In order to extend these approaches from matrix to tensor completion, we make use of extended notions on the concept of rank borrowed from linear algebra in conjunction with the high dimensional tensor theory. Tensor algebra is a mathematical framework that generalizes the concepts of vectors and matrices to higher dimensions. These tools are widely used to address problems of missing data in biomedical signal processing, computer vision, image processing, communication and seismic data processing. The High Order Singular Value Decomposition (HOSVD; Kolda, 2009) is used in many tensor rank-reduction algorithms. The HOSVD technique can be viewed as a generalization of the classical SVD for tensors. The computation of the HOSVD (De Lathauwer, 2000; Bergqvist, 2010) requires performing SVD as many times as the order of the tensor, over all the possible matrix representations of the tensor. Due to this fact the HOSVD can be expensive in terms of computational time. In tensor completion, the goal is to fill missing entries of a partially known tensor, under a low-rank condition. Many algorithms have been proposed in literature (Da Silva, 2013; Kreimer, GNGTS 2014 sessione 3.1